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Abstract 

We present an analytical study of electroosmotic flow in a Hele-Shaw configuration with non- 
uniform zeta potential distribution. Applying the lubrication approximation and assuming thin 
electric double layer, we obtain a pair of uncoupled Poisson equations for the pressure and depth- 
averaged stream function, and show that the inhomogeneous parts in these equations are gov¬ 
erned by gradients in zeta potential parallel and perpendicular to the applied electric field, re¬ 
spectively. We obtain a solution for the case of a disk-shaped region with uniform zeta potential 
and show that the flow field created is an exact dipole, even in the immediate vicinity of the disk. 
In addition, we study the inverse problem where the desired flow field is known and solve for 
the zeta potential distribution required in order to establish it. Finally, we demonstrate that such 
inverse problem solutions can be used to create directional flows confined within narrow regions, 
without physical walls. Such solutions are equivalent to flow within channels and we show that 
these can be assembled to create complex microfluidic networks, composed of intersecting chan¬ 
nels and turns, which are basic building blocks in microfluidic devices. 
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1. INTRODUCTION 


Electroosmotic flow (EOF) is the motion of a liquid due to interaction of an externally 
applied electric field with the net charge in the diffuse part of an electrical double layer. 
For solid surfaces with uniform zeta potential and chaimel geometry, EOF is character¬ 
ized by a uniform plug-like velocity profile. However, in practice, most surfaces are 
non-homogeneously charged to some extent, either due to manufacturing limitations or 
by intended design (see lUl and and references therein). 

EOF in capillaries with non-homogeneous zeta potential distribution has been studied 
thoroughly. Anderson and Idol |j3l found an exact solution to EOF through a capillary 
with axially varying zeta potential, assuming negligible inertia and a thin Debye layer. 
Herr et al. (T]] investigated, analytically and experimentally, EOF through a cylindrical 
capillary with an axial step change in zeta potential distribution. The authors measured 
experimentally the flow profile and showed good agreement with theoretical predictions. 
Ghosal in applied a lubrication approximation to investigate EOF in an infinitely long 
channel with slow axial variation in cross-section geometry and zeta potential. Ghosal 
showed that the flow rate through any section can be related to a uniform cylindrical 
capillary with an equivalent radius and zeta potential, which are determined solely by 
the geometry and surface charge distribution in the chaimel. 

Numerous works have also considered EOF between parallel plates. Ajdari |5ll3 was 
the first to analytically study and provide a closed form solution to the two dimensional 
problem of EOF between an undulating plate and a flat plate with a periodic surface 
charge distribution. Ajdari focused his analysis on the cross-section of the flow cell 
and demonstrated that the interaction of periodic deformations of the plate with peri¬ 
odic distribution of zeta potential gives rise to net flow generation between the plates, 
even though the plates are on average electro-neutral. Focusing on flow between parallel 
plates, Fong et al. |Zl provided an analytical solution for the three dimensional EOF field 
due to arbitrary distribution of zeta potential. In particular, by taking the limit of a small 
gap between parallel plates, the authors obtained a solution associated with a Hele-Shaw 
case, deducing that a localized defect in zeta potential distribution induces long-range 
flow perturbations. Ajdari |8l considered a three dimensional geometry, consisting of 
two plates with a spatially varying gap and periodic zeta potential in one direction. By 
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applying the lubrication theory, Ajdari formulated the Onsager matrix, which relates the 
applied pressure gradient and electric field to flow rate and electric current. Stroock et 
al. [91 performed an experimental study of EOF in flat shallow micro-chaimels having 
non-uniform zeta potential and compared the results with those of Ajdari (51 HI and Long 
et al. [Zl. Stroock et al. considered alternate positive and negative zeta potential pat¬ 
terns and demonstrated that various flow types such as multi-directional or circular flow 
can be generated, depending on whether the applied field is parallel or perpendicular to 
gradient of zeta potential. 

In this work, we aim to study the potential use of such non-uniform electroomotic 
flows as a mechanism to create complex microfluidic chaimel networks, not requiring 
solid walls. In sections II and III, we define the problem and derive a pair of uncou¬ 
pled Poisson equations for the pressure and depth-averaged stream function, governed 
by gradients in zeta potential parallel and perpendicular to the applied electric field, re¬ 
spectively. In section IV, we focus on axially symmetric zeta potential distributions, and 
show that the flow field created is an exact dipole in the immediate vicinity of a disk 
with uniform zeta potential. In section V, we show the effect of the orientation (relative 
to the electric field) of thin lines with uniform zeta potential on the pressure and vorticity 
fields. In section VI, we study the inverse problem where the desired flow field is known 
and solve for the zeta potential distribution required in order to establish it. We then 
demonstrate that such inverse problem solutions can be used to create flows of desired 
directionality confined within narrow regions, without physical walls. We show that 
these can be assembled to create complex microfluidic networks, and provide a specific 
example of Y-junction leading to a meandering charmel. 


II. PROBLEM DEFINITION 

We here denote dimensional variables by tildes and normalized variables without 
tildes. Figure 1 presents a schematic illustration of the geometry and the relevant phys¬ 
ical quantities. We consider creeping flow in a narrow gap between two parallel plates 
subjected to a uniform in-plane electrostatic field E\\. We employ a Cartesian coordinate 
system (x, y, z) whose x and y axes lie at the lower plane and 5 is perpendicular thereto. 
The gap between the lower and upper plate is h. Each has an arbitrary zeta potential dis- 
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tribution, respectively defined as C^(x, y) and which vary over a characteristic 

length scale I in the x — y plane. Hereafter, we adopt the || and ± subscripts to denote 
parallel and perpendicular vector components to the x — y plane, respectively 



FIG. 1. Schematic illustration of the problem, (a) We consider a configuration consisting of two parallel 
plates, separated by a gap h. Each of the plates is functionalized with arbitrary zeta potential distribution, 
(^{x,y), C^{x,y), on the lower and upper plates, respectively. The fluid is subject to a uniform electric 
field .E|j. (b) Averaging over the depth of the cell yields a 2D problem, in which denotes the average 
potential of the two plates, + C^)/2. 

The continuity equation for an incompressible fluid is 

V • ii = 0. (1) 


Assuming a thin electric double layer regime (see lElIll), the momentum equation in the 
bulk is given by 


P 



+ u ■ 'Vuj 


— Vp + 


( 2 ) 


where u = {it, v, w) is the velocity vector, i is time and p, ij and p denote the fluid density, 
dynamic viscosity and pressure, respectively. We account for the body forces acting on 
the double layer by using the Helmholtz-Smoluchowski slip boundary conditions, BTOll 


eC^Eu „ , eC^Eu 

'“II 1^=0 — ~ ) '“lllz=h ~ ~ ) '^Dz=0,h ~ 0; w) 

where itij = {u, v), u± = wz and e is the fluid permittivity, and are zeta potential 

distribution on the upper and lower plates, respectively. 


4 









III. GOVERNING EQUATIONS FOR EOF IN A HELE-SHAW CELL WITH NON-UNIFORM 
ZETA POTENTIAL 


Hereafter we denote characteristic scales of the system by an asterisk superscript. 
The typical magnitude of the flow velocity in the x — y plane, u*, is determined by the 
Helmholtz-Smoluchowski slip condition as u* = where (* is characteristic 

value of zeta potential and E* is characteristic externally applied electric field. The char¬ 
acteristic time scale is t* = l/E*, whereas the characteristic velocity in the z direction, w*, 
and the characteristic pressure, p*, remain to be determined from scaling arguments. 

We introduce the following normalized quantities, {x,y,z) = {x/l,y/l,z/h), t = i/i*, 
{u,v,w) = {u/u*,v/u*,w/w*), p = p/p*, C = C/C* arid E\\ = E\^/E*. We restrict our 
analysis to a shallow flow chamber, 

e = j < 1, (4) 


and a negligible inertia regime. The latter is characterized by a negligible reduced 
Reynolds number, eRe, defined as 


pu*h 

eRe = e^—:— 1. 

V 


(5) 


We assume that the Dukhin number, Du = as/cr^l, IHTIIT^ is small compared to e 


Du = <C e, (6) 

atl 

where ag/ab is the ratio of surface to bulk conductivities (see [121). As noted by Yariv, [I^ 
and Khair and Squires, [141 this ratio is a length scale (also known as the healing length) 
which determines the spatial variations in electric fields due to non-uniformity in surface 
conduction. High Dukhin numbers would result in spatial variations in electric field, as 
well as in variations in concentration which may lead to chemiosmotic flow corrections 
fl5] - [T7l . Under our assumption, Du e, the electric field and the bulk concentration can 
thus be considered to be uniform throughout the domain. 

From order of magnitude analysis of Q and (j^ we obtain, w* = eu* and p* = pu*/eH. 
We note that as the flow is driven by the wall boundary conditions, the characteristic 
pressure is independent of the viscosity, p* = —eC*E*/eR. Expanding the velocity field 
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and the pressure in powers of e, llTSlI and considering the leading order, 0(e°), the nor¬ 
malized continuity and momentum equations Q and (j^ take the form 


du dv dw 
dx dy dz 


(7a) 


V||P=^ + 0(ei?e,e2), (7b) 

^ = 0{e^Re,e^). (7c) 

Integrating ( [7b] > twice with respect to z while making use of the boundary conditions Q, 
we obtain an expression for the in-plane velocity field 


u 


(2 -i)V||p+^ £;ii + c"'£;ii. 


( 8 ) 


where Vy = {d/dx, d/dy) is the two dimensional gradient. Similarly, the perpendicular 
velocity can be expressed as 


Ux = z{l - z)E\\ ■ [^V||C^ + (^ - 1)V||CT 


(9) 


Defining the mean in-plane velocity, as (u\<^ = u\\dz, and making use of (j^, (7a) and 
^ yields 

V||-(u|i> = 0, (10) 


and 

(“ll> = -^V||p+{0-B||. (11) 

where {() is an arithmetic mean value of the zeta potential on the walls, {() = + (^)/2. 

Applying the two dimensional divergence to ( [TT] ), and using ( [lO] ), we obtain an equation 
in terms of the pressure only, 

Vjp = 12E^^ ■ V|| (C). (12) 

Similarly, applying the normal component of the curl operator to ( [TT] ), leads to an equa¬ 
tion for the stream function. 


= {E^^ X V|| (0) ■ 2, (13) 

where ^jJ{x, y) is the averaged stream function related to the velocity field through ) = 

{di/j/dy, —d^j/dx). 
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We note that (12) and (13) admit an associated gauge freedom in the choice of zeta po¬ 
tential, which does not affect the resulting pressure or the flow field. Specifically, for the 
case of an electric field acting in the x direction, it follows that adding an arbitrary func¬ 
tion Cipiy)^ to the zeta potential, {({x, y)) —)■ {({x, y)) + Ci>{y), will not modify the resulting 
pressure in ( [l^ . Similarly, the stream function in ( [l^ is indifferent to the transformation 
{C,{x,y)) {C,{x,y)) + Cp(2:)- 

Eqs. ( [T^ and (13) are an uncoupled set of Poisson equations for the pressure and the 
stream function, and extend the Hele-Shaw equation IIT9l to include non-uniform EOF. 
Notably, ( [T2] > includes a source term that depends on gradients of zeta potential which 
are parallel to the applied electric field, while (131 relates the vorticity, u = — V|'0, to 
changes of zeta potential in a direction normal to the applied electric field. In particular, 
in a region where the non-homogeneous term in ( [T^ vanishes, the resulting flow field 
would be irrotational and a velocity potential function could be defined. 

In addition to direct calculation of the flow field arising from a certain zeta distribu¬ 
tion, we can also use ( [T^ to obtain the required zeta potential distribution for achieving 
a desired flow field. Without loss of generality, we may assume that the uniform electric 


field is directed along the x axis, E\\ = Ex. Solving (13) in terms of a desired stream 
function we obtain the zeta potential distribution 


{C{.x,y)) = ^ I ^l^{x,y)dy + Cpix), 


(14) 


where Cp{x) is an arbitrary function of the x argument to be determined from the bound¬ 
ary conditions. 


IV. AXIALLY SYMMETRIC ZETA POTENTIAL DISTRIBUTION 

Consider the case where the zeta potential distribution acquires a non-zero value, Co/ 
in the inner region (ro < r) and vanishes in outer region, (ro > r). In the following, we 
adopt the superscripts in and out to distinguish between physical quantities in each one 
of the two regions. The corresponding zeta potential distribution is given by 

(C(r)) = Coi^(ro-r), (15) 

where r is the distance as measured from the center of our coordinate system which coin¬ 
cides with the center of the disk, H stands for the Heaviside step function, tq is the radius 
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of the disk, which here represents a characteristic length scale /, and Co is the constant 
value of the zeta potential in the disk. 

Clearly, the abrupt change in zeta potential value described by ( [T5] > locally violates 
the lubrication approximation, which assumes negligible gradients in the x — y plane. 
Nevertheless, the resulting deviations are expected to be limited to a narrow region, of 
order e, similar to deviations near the step changes of a cross section in a long and narrow 
channel 


For = {E, 0), (12) takes the following form in polar coordinates 


1 d f dp\ 1 d‘^p 


“7pm 7m I + 

r or 


— 12E cos{9) 


d{0 


= 0 , 


(16) 


dr J dO'^ ^ dr 

where 6 is the azimuthal angle in the plane and x = r cos{9). The last term in ( [1^ suggests 
a solution of the form 

p{r, 9) = a{r) cos{9), (17) 

where a(r) is yet to be determined. Combining ([T^-(|l7|) yields 


Id/ da{r] 
rdr 


^ + 12ECo6{ro - r) = 0, 


(18) 


dr I r^ 

where 5 is the Dirac delta function. Solving ( [l^ separately in the iimer and outer regions, 
and requiring regularity yields the following expression for the pressure 


p(r, 9) = 



_|_ cos(6*) 


(19) 


E'^r cos{9)^ 


where 6"“*, and the dipole strength are coefficients to be determined from bound¬ 
ary conditions. Utilizing (11), (15) and ( [T^ provides the associated closed-form expres¬ 
sion for the flow field 


(^ll> = 



Lout 


cos{9)r + 


^out 


+ b 


,out 


sin(9)0 


— (-E^ + 12UCo) cos(9)r + (E^ - 12UCo) sm(9)0 


ro < r 


r < tq. 


( 20 ) 


In the iimer region, the pressure changes linearly in the direction of the electric field, 
which results in a uniform velocity vector field, explicitly given by 




X. 


( 21 ) 
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The coefficients are determined by demanding continuity of the pressure and radial ve¬ 
locity component at r = tq, leading to 

= QECorl' = 6 ^( 0 - ( 22 ) 


The dipole strength 0 °“* represents a contribution to the pressure as a result of the zeta 
potential in the disk, while the constant represents the magnitude of a uniform flow 
field in the x direction which stems from an externally applied pressure gradient in that 
direction, = Ap/Ax (e.g. tor the case of constant pressure at r —)■ 00 , = 0). 

is then readily obtained from (22). The value of 6 *” — represents discontinuity in 
the 0 component of the depth-averaged velocity on the surface r = tq, induced by the 


vorticity sources specified by (13). From (20) we find that the corresponding jump in the 
6 component of the velocity (i.e. difference between its outer and inner values) is —E(o. 


This result can be also obtained by using the relation {ue) = —d'llj/dr and integrating (13) 
along infinitesimal region around tq. 

While it is expected from multipole expansion theorem II 2 TII that far from the non¬ 
uniformity the flow field will decay to a dipole IZl, we note that in this case the flow 
exhibits an exact dipole behavior even in immediate vicinity of the disk. 

Figure (2a) presents the streamlines and pressure distribution map of a single elec- 
troosmotic dipole for the case of vanishing pressure far from the disk (6°“^ = 0). As 
expected, the boundary conditions in the inner region dictate that the flow is from low 
pressure to high pressure while the flow in the outer region is from high pressure to low 
pressure. Interestingly, taking advantage of the uniform flow field in the iimer region, we 
can superpose it with another depth-averaged flow of equal magnitude and opposite di¬ 
rection, to achieve zero net flow in the inner region, i.e. = 0. This can be achieved 

by applying a mean pressure gradient Ap/Ax = = 6E(q, or by adding a bias value 

of —Co/2 to the zeta potential everywhere in the domain. Figure (2b) presents the stream¬ 
lines and pressure distribution for this case. The corresponding flow field around the 
disk coincides then with the well-known solution of potential flow past an infinitely long 
cylinder. 

Employing the linearity of the governing equations enables to superpose the basic 
dipole solution, ( [T^ , and construct axially symmetric solutions where the zeta potential, 
Co(r), is a function of the radial coordinate. For the case in which Co(^) vanishes outside a 
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FIG. 2. The pressure distribution (colormap) and streamlines (white lines) due to a uniform electric field in 
the X direction, applied to a Hele-Shaw configuration having a non-zero zeta potential in a disk of radius tq. 
(a) The case of no incident flow. The depth-averaged flow field described by ^0\ , results in a uniform flow in 
the inner region and dipole flow in the outer region. Note that the points of extremum pressure and vorticity 
correspond to locations where the source terms in ( [l^ and ( [l3| attain their minimal or maximal values, (b) 
The case where an external pressure gradient (counter-flow) is applied such that the velocity inside the disk 
vanishes. The result coincides with that of potential flow around a cylinder and is characterized by open 
streamlines which do not enter the disk region. Both solutions are normalized such that the pressure varies 
between ±1. 


disk of radius vq, we obtain the corresponding pressure distribution, ptot, 


Plot {r, 9) = 


12Ecos(0) f^ , 


12E cos{9) 


I r ^ /ro 

- I CoiroYodro + ^ U YroYod^o 

0 ^ \r 


Vq < r 
cos(6*) r < ro, 


(23) 


where the solution in the outer region corresponds again to an exact dipole, having a 

ro 

dipole strength 12£' j Co(?"o)’"o'^^o- Figure (3) presents an analytical solution for superpo- 

0 

sition of two concentric disks of radii ro/2 and tq with zeta potentials Co arid — Co/4, re¬ 
spectively. For this special case, the pressure vanishes outside the larger disk, and results 
in flow that is confined to the boundaries of the outer disk, whereas a uniform velocity 
field is obtained in the small disk, as expected. 

One can also readily apply a combination of disk-shaped regions with constant zeta 
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FIG. 3. Analytical solution due to superposition of two disks, described by Eqs and ( [^ , showing the 
superposition of two disks of radii ro/2 and ro with zeta potentials, Co (ind —Co/4, respectively. Colormaps 
describe the normalized pressure distribution and the white lines show streamlines. The resultant dipole 
strength vanishes outside the larger disk, giving rise to a pressure and flow which are confined within the 
outer disk and uniform flow in the inner disk. 



potential as a tool to achieve a required flow field in a Hele-Shaw cell with a uniform 
electric field. Specifically, as an illustration of this approach, we apply the panel method, 
commonly used in aerodynamics Il23 . to design EOF around a symmetric NACA 0015 
airfoil profile. For concreteness, we choose the direction of the electric field in our exam¬ 
ple along the x axis, and specify a known pressure gradient, Ap/ Ax, along this direction. 
We place a set of 24 disks, each having a uniform (but potentially different) zeta poten¬ 
tial, also along this axis. We then utilize the closed form dipole solution generated by 
each disk to calculate the velocity vector at 24 points along the airfoil curve. Demanding 
no-penetration at each of the points, results in a set of linear algebraic equations for the 
associated ro and Co values of the disks. 

Figure (4a) shows the intensity of each of the disks as well as their location along the 
X axis. Figure (4b) presents the streamlines (solid lines) obtained from the analytical 
solution. To validate the results, we performed a three dimensional direct numerical 
simulation using Comsol Multiphysics 4.4, in which the Stokes flow equations are solved, 
coupled to Helmholtz Smoluchowski boundary conditions. Our grid consisted of 120,504 
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Numerical 
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FIG. 4. Superposition of disks with uniform zeta potential enables creation of complex flow patterns, satis¬ 
fying no-penetration on desired convex surfaces, (a) 24 circular disks are distributed along the xaxis . The 
non-dimensional dipole strength of each of the disks is set such that a no-penetration boundary condition is 
obtained over a NACA 0015 airfoil, (b) Comparison of analytical (solid lines) and numerical results (dashed 
lines) of the flow field obtained from the distribution. Analytical results were obtained from 2D Hele-Shaw 
model, whereas the numerical results are depth-averaged from a direct three dimensional numerical solution 
using Comsol Multiphysics 4.4. 

2 nd order unstructured prism elements, and all solutions converged by at least seven 
orders of magnitude from a zero velocity initial condition. We used equally sized disks 
having zeta potential values as determined by the analytical solution. The velocity field 
is depth-averaged, and streamlines of this numerical solution are also presented in the 
figure (dashed lines), showing good agreement with the analytical result. 


V. UNIFORM ZETA-POTENTIAL ALONG A THIN FINITE LINE 


Eqs. ( [l^ and (131 are both Poisson-type equations for the pressure and vorticity, in 
which gradients in zeta potential serve as source terms. The respective dot-product and 
curl operators in these source terms suggest that the solutions would depend on the ge¬ 
ometry of surface patterning. To further highlight the roles of these source terms, let us 
consider two simple cases of EOF arising from a thin finite line with uniform zeta poten¬ 
tial which is either aligned with the electric field or perpendicular to it. In both cases, we 
model the lines using the linearity of the governing equations, and sum up the pressure 
contributions of dipoles having a strength per unit length 


dp{x,y) = 


a°^*{x — Xo{l)) 


[x - Xo{l))‘^ + {y - yoil)y 


dl, 


(24) 
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where the dipole position is parametrized by I, and dl is an infinitesimal interval along 
the curve 1. We then use ( |TT| > to derive the fluid velocity profile outside the dipole distri¬ 
bution. 

For the case of a line of dipoles direcfed along the electric field, the pressure distribu¬ 
tion and the mean velocity are given by 


pix, y) = 0°“* In — 
r, 


r_ 


(25a) 


(^ll> = 


-.out 


12 


-)■ 

r+J 


(25b) 


where f± are unit vectors pointing from the observation point (x, y) to the edge points 
(x_,|/o) and (x+,|/o) , and r± = \J{x — x±)‘^ + (y — yof are the corresponding distances 
between these points. For a case where the line coimects the points (xq, y-) and (xq, y+) / 
and is perpendicular to electric field, the pressure distribution and the mean velocity are 
given by 


p{x,y) = a 


out 



-1 ( y-y+ 

X — Xq 


(26a) 

(26b) 


where 6± are unit vectors which perpendicular to f±. 


The solutions (25a I and (25b) can be interpreted as a sink and a source, located at the 
edges of the uniform zeta potential line. Figure (5a) presents the resulting pressure and 
velocity field for this case. Clearly, the apparent violation of continuity, V||- (u||) = 0, at 
the two points is an artifact of the fact that ( 25a[ ) and ( 25b[ ) represent the flow outside the 
array of dipoles, and additional flow that coimects (x_, yo) and (x+, yo) exists within this 
array. Note that in this case all zeta potential gradients are in the direction of the electric 
field and thus, consistent with the streamlines equation ( [T^ , the flow field ( 25b[ ) is free 
of vorficify. The solutions (26a) and (26b) are effectively equivalent to a pair of equal 
strength and opposite sign vortices, centered at the edges of the dipole line. Figure (5b) 
presents the exact analytical solution for the pressure and velocity, as given by ( 26a[ ) and 
(26b). Consistent with ( [l^ , high pressure gradients are formed befween the two sides of 
the line, where the zeta potential changes in a direction parallel to electric field. At the 


line's edges the pressure difference results in circulatory flow, again consistenf wifh (13) 
which predicts vorticity where the zeta potential changes in a direction normal to electric 
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field (in this case, the uniform zeta line terminating into the zeta-free surface). The circu¬ 
lation, defined as the line integral, T± = ^ (un) ■ dZ, along a closed path around each edge 
point coincides (up to a factor of 7r/6 ) with individual dipole strength r± = ±TTa°^^/6. 


Clearly, r+ -|- r_ =0, i.e., the total circulation remains zero. More generally, from (13) 
it follows that the total vorticity generated by non-homogeneity of an arbitrary shape 
which hosts a constant zeta potential must vanish: the component of the zeta potential 
gradient perpendicular to the electric field is {() sin(6*), and thus its line integral along the 
closed curve must vanish. 



FIG. 5. Analytical solution showing the flow field due to a line of uniform zeta potential aligned (a) along 
the direction of the electric field and (b) perpendicular to the electric field. Each solution is obtained from 
superposition of dipoles, according to ( [^ . Colormaps describe the normalized pressure distribution and 
the white lines show streamlines. In configuration (a) gradients in zeta potential are all in the direction 


of the electric field, resulting in source/sink terms at the edges of the line for the pressure equation (12). 
The streamline equation ( (l3] >, on the other hand, remains homogenous and no vorticity is formed. In (b), 
the lines edges form a zeta potential gradient perpendicular to the electric field, resulting in vorticity and 
localized circulation around each for the edges. In addition, a pressure difference forms across the line, 
reflecting the fact the zeta potential changes in a direction parallel to the direction of the electric field. 
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VI. ZETA POTENTIAL DISTRIBUTION DEFINED BY A REQUIRED FLOW FIELD AND 
ITS APPLICATION TO CREATE COMPLEX MICROFLUIDIC NETWORKS 


In sections IV and V the zeta potential distribution geometry was predefined and then 
solved to obtain the flow field. More complex flow fields were obtained by superposing 
solutions. Here we focus on obtaining the zeta potential distribution required in order 
to establish a desired flow field. This approach enables to generate depth-averaged flow 
fields which caimot be created using external pressure actuation at the boundaries. This 


is done by setting a desired stream function, and using (14) to calculate the required 
zeta potential distribution. Without loss of generality, we here assume that the uniform 
electric field is directed along the x axis, E\\ = Ex. 

We illustrate this approach by considering flow in a square domain [0,1] x [0,1], where 
X = 0, 1 and y = 0 correspond to impermeable walls, and y = 1 corresponds to a surface 
of constant pressure (see figure (6a)). An example of such a flow is given by the stream 
function, 

= l/sin(7rx), (27) 


describing a U-shaped flow (see solid streamlines in figure (6a)) in which the flow enters 
and exits the same y = 1 face. Utilizing ( [14| ) and ( [27[ ) yields the corresponding expression 
for the zeta potential 

(C(a^, y)) = -^2/^sin(7rx) + Cp(a:). (28) 

Defining the normal and tangential unit vectors to the domain boundaries as n and i, re¬ 
spectively, the no-penetration boundary condition, ■ f = (un) ■ n = 0 is automatically 
satisfied by the stream function ( [27| ) on x = 0, 1 and y = 0. To satisfy the constant pres¬ 
sure condition, Vyp ■ f = 0, on the face y =1, we utilize the gauge freedom in the choice 
of the function Cp(x). Substituting this condition into ( [TT| results in (u\\) ■ x = {Q \y=i E. 
Expressing the velocity in terms of the stream function, and substituting the expression 
for zeta potential ([^ leads to 


Cp(^) 


1 

E 



sin(7rx). 


(29) 


Consequently, the corresponding zeta potential distribution and the pressure are 





sin(7rx) -|- 2] sin(7rx), 


(30a) 
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p{x, y) = Gvr {y^ — l) cos(7ra;), 


(30b) 


and are presented in figures (6b) and (6c), respectively. 


To validate the results, we compare the desired flow field (271 with the flow field ob¬ 
tained by numerically solving a three dimensional case in which the resulting zeta poten¬ 
tial ( 30a[ > is used in the slip boundary conditions on both the lower and upper walls. To 
solve the corresponding Stokes equations, we used the creeping flow package in Comsol 
Multiphysics 4.4. Streamlines of the depth-averaged velocity field are presented in figure 
(6a) (dashed lines), showing good agreement with the desired flow. 


(a) Streamlines (b) C -potential (e) Pressure 



Numerical E 




FIG. 6. Demonstration and numerical validation of the inverse problem, in which the zeta potential distri¬ 
bution is calculated from a desired flow field, (a) Comparison of the desired streamlines (solid lines) with 
the depth-averaged streamlines obtained from a three dimensional numerical simulation (dashed lines), (b) 


The analytical zeta potential distribution, \30a\, obtained in order to generate the desired flow field, (c) The 


resulting pressure distribution {30b). All calculations are performed using E = 1. 


We now turn to illustrate the use of ( [14] > to engineer complex flow fields of the type 
encountered in microfluidic devices. These flows are typically confined to narrow chan¬ 
nels composed of turns, junctions and straight segments. Importantly, in our approach 
there are no real walls that confine the liquid, and yet using appropriate zeta potential, 
the fluid achieves desired directionality within specific regions of interest and its velocity 
strongly decays outside them. We show that under certain conditions the corresponding 
zeta potential distribution necessary to generate complex flow fields of such type admits 
modularity, i.e. can be represented as a sum of simple and basic distributions which allow 
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sufficiently smooth matching. 

First we wish to determine the zeta potential necessary to generate a flow field along a 
curve y = f{x), under a homogeneous electric field. We set a stream function of the form 

^l^{x,y) = F{y-f{x))^F{x), (31) 


which ensures that xjj has a constant value along y = f{x) for any choice of F. The variable 
X = y — f{x) parametrizes different streamlines according to their displacement along the 
y axis. By definition of the stream function, and using ( [3T] >, the velocity components can 
be expressed as 

= =/'WF'(x). (32) 

To determine the zeta potential necessary to generate the flow field ( [3^ , we substitute 
the stream function (j3^ into the relation ( [T4] > which yields the following closed form 
expression for the necessary zeta potential 


= g (F{x} + f[xfF'{x} - F(x)f"(x))+Ux), 


(33) 


where (p(x) is an arbitrary function. Hereafter we focus on a case where (p(x) = 0. Uti¬ 


lizing (32) and polar representation of the flow field components, (u) = (uy) cos( 6 ') and 
(v) = (uii) sin( 6 '), (33) is rewritten as 




(34) 


This expression allows to directly obtain the zeta potential required in order to pro¬ 


duce the desired flow, described by the stream function (31). In regions where the flow 
direction changes significantly the second term in (j3^ dominates over the first term, 
and vice versa in regions where the flow field is nearly unidirectional. In particular, 
in regions with uniform flow the zeta potential admits a simple expression given by 


(u) /Ecos^{6) = (u\^ /Ecos{6). Also, (34) shows that the zeta potential required to gen¬ 
erate flow in a direction perpendicular to electric field (i.e. (m) = 0 ), diverges and is thus 
non-physical. 

A convenient function to define a flow field in which the velocity is confined to a 
narrow region around the streamline y = f{x) is 


ij^x, y; 9in, Oout) = tanh (7 {y - f{x; 9in, 9out)) ■ 


(35) 
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This function exhibits significant gradients in a region of width I/7 around the stream¬ 
line. Outside of this region it quickly reaches a negligible gradient, and thus negligible 
flow. 

For creation of stream function modules, which are useful for the construction of mi- 
crochaimel segments, it is convenient to choose the function / as 


Bini dout) tan 


X 


tan {9out)- 


X 


(36) 


(1 + exp (/5x))" ^ (1 + exp (-/3a;))" 

This function describes a smooth approximation of a bi-ltnear curve, which approaches 
y = tan {9out)x as x tends to -|-cxd and io y = tan {9in)x as x tends to —cxd . n and /3 are 
positive parameters which set the exact shape of the curve near the bending point ( 0 , 0 ). 

Figure (7a) presents the contour map of the stream function, - 0 ( 0 ;, y; — 7 r/ 4 , 0). The 
streamlines describe an incoming fluid flowing at an angle — vr/4, a turn region around 
the origin and outgoing flow along the x axis (zero angle). The region which supports 
a significant flow velocity is enclosed between the dashed lines, and its width is of or¬ 
der 1 / 7 . The magnitude of the flow velocity on the dashed lines drops to 1% of its peak 
value. Figure (7b) shows several streamlines (blue arrows), together with the zeta poten¬ 
tial described by ( [M) >, necessary to generate this flow field. The second derivative of the 
function ^(a:, y; — 7 r/ 4 , 0) acquires significant values around the line x = 0, which results 


in a significant contribution of the —F{x)f"{x) = —'ipix, y)f"{x) term in (33). Since f"{x) 
is independent of y, and since the stream function values are constant far from the chan¬ 
nel region, this results in zeta potential values of opposite signs at x = 0 , which extend 
uniformly toward y = ±cxd. Far from the line x = 0, where the flow has a nearly uniform 


directionality, the first term in (34) dominates and the resulting zeta potential distribution 
is uniform along the incoming and outgoing flow directions. 

More complex flows can be created by joining of multiple stream function modules. 
Figure ( 8 a) shows the stream function, 1/2 ('0(x, y; —vr/4, 0) — 'ijj{x, y; -l-vr/4, 0)), obtained 
by superposition of two stream functions described by ( [3^ . In such superpositions, the 
resulting flux is thus the sum of individual fluxes. This allows, for example, construct¬ 
ing junctions where multiple incoming streams are combined into a single one. In the 
example shown in figure ( 8 a), the resulting flow consists of incoming streamlines along 
the angles — 7 r /4 and 7 r /4 towards the origin, and a combined outgoing flow along the 
X axis. We note that the total flux in the chaimels can be controlled by modifying the 
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FIG. 7. Analytical example ofzeta potential calculation for a bent channel segment, (a) Colormap of the 


stream function 'f{x,y,—7r/A,0) given by {351 which describes incoming flow at an angle —TrjA and 
horizontal outgoing flow. Here /3 = 7 = 2 , re = 4, = 1 . The stream function is constructed of 

two regions having negligible gradients (and thus negligible flow), which are connected through a steep 
gradient region of width I/ 7 . The dashed lines correspond to streamlines on which the velocity drops to 1% 
of its maximum value, thus defining the outer boundaries of a channel, (b) Colormap of the zeta potential, 
necessary to generate the desired flow field (blue streamlines), obtained by inserting the corresponding 


stream function into {33 iThe pattern shows two regions of opposite values ofzeta potential which extend 
along the line x = 0, towards positive and negative directions of the y axis. This contribution stems from the 
term —F{x)f"{x) in ( |3^ , and occurs in the region where f(x) admits its highest second derivative. Two 
other regions which host significant values ofzeta potential correspond to zones with nearly unidirectional 
flow, represented by the first term in 


value of whereas the relative flux of each of the streams can be dictated by 

modifying the value for fixed values of and Figure ( 8 c) shows the zeta 

potential necessary to generate this flow field, obtained by substituting the correspond¬ 
ing stream function into ( [34] >. Note that following the same steps one can consider a linear 
combination of N basic functions fj{x, y] 9n, 0) which generates a converging flow from 9n 
incoming directions and a combined flow along a single outgoing direction. Figure ( 8 b) 
presents streamlines obtained from matching the stream functions — Xi,y; —7r/4, 0) 
and fj{x — 0 : 2 , 2 /;+ 7 r/ 4 , 0), (with X 2 = —xi = 3.5), along the common boundary x = 0. 
Matching allows the construction of extended chaimels in which the flux in one module 
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Superposition of modules 



-4 0 4 ™ -4 0 4 ™ 

Matching of modules 



FIG. 8. Analytical examples showing superposition and matching of two basic modules, (a) Superposition 
assembly of the functions l/2V’(x, y; —7r/4,0) and forming a Y-junction in which 

incoming flows are directed at an angle ±7r/4, and combine into a single outgoing flow along the x axis. 
The value determines the total flow rate in the channel, whereas modifying the value of'f^^'i 

allows to control the fraction carried by each of the two incoming flows, (b) Stream function of two modules: 

—vr/4, 0)forx < Oand 'f{x—X 2 , y; +7r/4, 0)forx > 0 (with X 2 = —xi = 3.5), matched along 
the line x = 0, forming U-shaped streamlines. (c,d) Colormaps of the zeta potential distributions necessary 
to generate the desired flow field (blue arrows), obtained by substituting the stream functions corresponding 


to (a,b) into (33 In all subfigures the magnitude of flow velocity on the dashed lines drops to 1% of its peak 
value and the parameters which enter ^3, ([^ and ^3 are given by (3 = y = 2, n = 4, E = 1. 


is preserved in the next. In this assembly, the resulting flow consists of an incoming 
stream along a straight line with slope — tarL(27r/5), a short horizontal flow parallel to 
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the X axis, and outgoing flow along a straight line with a slope tan (27r/5). To join two 
such elements, the values of C arid its derivative along the common boundary should 
ideally be matched. ( [M] > shows that the zeta potential of streamlines directed along the 
electric field (zero slope), depends solely on the velocity value and angle and thus can 
be used for outgoing flow from one domain and incoming flow to another. The resulting 



FIG. 9. Analytical example demonstrating the use of matching and superposition of basic modules for 
creation of a microfluidic network, (a) and (b) show the desired streamlines pattern and the zeta po¬ 
tential distribution necessary to generate the corresponding flow field, respectively. The microfluidic 
network is composed of a Y-junction, created by superposition assembly of l/2?/)(x, y;+7r/4,0) and 
—l/2ip{x, y; —vr/4, 0), followed by a matching assembly of — Xi,y — fi; —7r/4,0) and 'ijj{x — Xi,y — 
fi] +7r/4,0) where the index i enumerates the bending points (xj, f) of different basic modules, and their 
matching boundary lines, x = (x* + Xi+i)/2. Detailed structure of each one of the modules enclosed within 
the regions (i), (u), is described in figure (8) while region (in) is described by a module similar to the one 
shown figure (7). In both colormaps the magnitude of flow velocity on the dashed lines drops to 1% of its 
peak value and the parameters which enter ([3^ and ([3^ are given by f3 = y = 2, n = 4, E = 1. 
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zeta potential for this type of matching is presented in figure (8d). While the values of 
the stream function, and zeta potential match precisely along the common boundary, 
X = 0, their corresponding derivatives experience a discontinuity The discontinuity in 
the stream function derivative, A , along the common boundary is given by 


^ ^ 1 + (1 — n/3x) exp (dx) 

(l + exp(/?x)r+' 


(37) 


and could be made arbitrarily small by extending the length of the domain or by proper 
choice of the corresponding parameters. For typical values in our example, A does not 
increase above 10“®, which would be negligible in typical microfluidic systems. 

Figure (9) illustrates the use of superposition and matching for creation for a mi¬ 
cromixer geometry, in which two streams merge in a Y-junction, and continue to a serpen¬ 
tine chaimel consisting of six 47r/5 turns. Figure (9a) shows the stream function, which is 
an assembly of basic flow fields matched together along common boundaries. The mag¬ 
nitude of flow velocity on the dashed lines drops to 1% of its peak value. Figure (9b) 
shows the necessary zeta potential distribution, obtained from ([M)>. 


VII. SUMMARY AND CONCLUSIONS 

In this work we studied EOF in a Flele-Shaw configuration, with non-uniform zeta 
potential distribution. We demonstrated the use of zeta potential distributions for pat¬ 
terning of complex flow fields, including microchaimel networks. 

The flow fields we obtained in this work are depth-averaged fields. Care should be 
taken when analyzing the flow of particles in such fields; for sufficiently high diffusivity 
of the particle (i.e. low Peclet number), an ensemble of particles will move with the 
average velocity field. However, at high Peclet, individual particles may follow particular 
streamlines and the ensemble will exhibit a behavior substantially different than that of 
the average. 

In sections IV-V we demonstrated the construction of complex flow fields by super¬ 
position of basic solutions. We note that ( [TT] >, ( [l^ and ( [l^ are conformal invariant. This 
follows directly from the fact that these equations can be rewritten so that each term 
will contain equal number of gradient operators (see 112311 1. This opens a door to em¬ 
ploy conformal mapping machinery, to calculate resultant pressure and flow field due to 
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non-uniform zeta potential patterning in the presence of complex geometries with real 


walls. It is also worth noting that the governing equations (12) and (13) tor the pres¬ 
sure and depth-averaged stream function can be extended to include the effect of slow 
varying height between the plates. The modified governing equation will be an EOF 
driven Reynolds equation in which both the viscous terms and the EOF terms depend on 
the spatially varying gap between the plates. While the governing equation will remain 
linear, and thus could potentially be resolved analytically, it is beyond the scope of this 
work. 

In practice, surface patterning could be obtained in several ways, including chemical 
modification of the surface, or using an array of surface electrodes. The latter adds com¬ 
plexity, but has the potential advantage of allowing surface potentials which are much 
higher, as well as dynamic modification of the zeta potential distribution. Further study 
would be required to characterize the effect of zeta potential discretization on the result¬ 
ing flow fields. Flow patterning by EOF, particularly by dynamic modification of the 
surface potential, may serve a powerful tool for manipulating fluids without mechanical 
components in a variety of microfluidic applications. 
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